# 2022.3.23
# re-analyze smart-seq log2(TPM+1) data with Seurat
# from 2020 Zhang Cell paper 

# edit in 2023/5/30

library(data.table)
library(Seurat)
library(tidyverse)
library(tidyseurat)

# load smart-seq TPM data ---------
data <- fread('CRC-I/data/GSE146771_CRC.Leukocyte.Smart-seq2.TPM.txt.gz')

data <- column_to_rownames(data, 'V1')

II_smart <- c('P0104','P0411','P1212')
IT_smart <- c('P0305','P0720','P0825')

data <- select(data, contains(c(II_smart, IT_smart)))

# load metadata
smrt_meta <- fread("CRC-I/data/GSE146771_CRC.Leukocyte.Smart-seq2.Metadata.txt.gz")

## recover counts matrix ------
# from nUMI and genelength
cell_umi <- tibble(CellName = colnames(data))

sorted_umi <- cell_umi |>
  left_join(smrt_meta) |>
  select(CellName, raw.nUMI, filter.nUMI)

## create seurat ------------
sobj.smt <- data |>
  CreateSeuratObject(min.cells = 3,
                     min.features = 200,
                     names.field = 3)

# assign I232T genotype in smrt_meta
sobj.smt <- sobj.smt |>
  mutate(CellName = .cell,
         genotype = case_when(
           orig.ident %in% II_smart ~ 'II',
           orig.ident %in% IT_smart ~ 'IT',
           .default = 'NA'
         )) |>
  left_join(smrt_meta)

sobj.smt |>
  FindMarkers(ident.1 = 'hM12_TAM-C1QC', ident.2 = 'hM13_TAM-SPP1', group.by = 'Sub_Cluster', features = 'C1QB')

FetchData(sobj.smt, 'TPSAB1') |>
  summarise(max(TPSAB1))

sobj.smt |>
  ggplot(aes(nFeature_RNA, filter.nGene)) + geom_point()

sobj.smt$mitoRatio <- PercentageFeatureSet(object = sobj.smt,
                                           pattern = "^MT-")

sobj.smt |>
  VlnPlot('mitoRatio')

# save genotyped data
write_rds(sobj.smt, 'CRC-I/data/smart_tpm.rds')

# load 10x TPM data ------------
data <- fread("CRC-I/data/GSE146771_CRC.Leukocyte.10x.TPM.txt.gz")

data[1:5, 1:5]

data <- column_to_rownames(data, 'V1')

II_10x <- c('P0410','P0323','P0408','P1026','P0104')
IT_10x <- c('P0123','P0202','P0613','P1025','P0305')

# the cell naming: P_N_P0410_01391, first letter is always 'P', second letter is N/P/T representing tissue source.
# here we only need cells in tumor
data <- dplyr::select(data, contains(c(II_10x, IT_10x)))

# load metadata
tenx_meta <- read_delim("CRC-I/data/GSE146771_CRC.Leukocyte.10x.Metadata.txt.gz")

## recover counts matrix ------
cell_umi <- tibble(CellName = colnames(data))

sorted_umi <- cell_umi |>
  left_join(tenx_meta) |>
  dplyr::select(CellName, raw.nUMI, filter.nUMI)

data[1:5,1:5] |>
  map2(sorted_umi$filter.nUMI[1:5], \(x,y)x*y/10000)

sobj <- CreateSeuratObject(data,
                           names.field = 3,
                           min.cells = 3,
                           min.features = 200)

# assign I232T genotype in tenx_meta
sobj <- sobj |>
  mutate(CellName = .cell,
         genotype = case_when(
           orig.ident %in% II_10x ~ 'II',
           orig.ident %in% IT_10x ~ 'IT',
           .default = 'NA'
         )) |>
  left_join(tenx_meta)

sobj$mitoRatio <- PercentageFeatureSet(object = sobj, pattern = "^MT-")

sobj |>
  VlnPlot('mitoRatio')

sobj <- sobj |>
  FindVariableFeatures() |>
  ScaleData() |>
  RunPCA() |>
  RunHarmony('orig.ident') |>
  RunUMAP(reduction = "harmony", dims = 1:20) |>
  FindNeighbors(reduction = "harmony", dims = 1:20) |>
  FindClusters()

VlnPlot(sobj, 'ACTB', pt.size = 0, group.by = 'Global_Cluster')

sobj <- mark_cell_type_singler(sobj, ref = hpca)
DimPlot(sobj, group.by = 'Global_Cluster')
DimPlot(sobj, group.by = 'singler_label')

# save 10x seurat file
write_rds(sobj, 'CRC-I/data/zy2020_tumor10x.rds')
